Distance preconditioning for lattice Dirac operators 
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We propose a preconditioning of the Dirac operator based on the factorisation of a predefined 
function related to the decay of the propagator with the distance. We show that it can improve the 
accuracy of correlators involving heavy quarks at large distances and accelerate the computation of 
light quark propagators. 



I. INTRODUCTION 



A key ingredient of lattice QCD simulations is the in- 
version of the Dirac operator which enters the generation 
of unquenched gauge field configurations and the compu- 
tation of hadronic observables. One needs to solve nu- 
merically a linear system of the form 



J2(D[U}+M) x , y S(y)=r)(x) 



(1) 



where D[U] is the chosen discretization of the massless 
interacting Dirac operator, M is the quark mass in lattice 
units, Tj(x) is a source vector that is different from zero 
on a single time-slice (that without any loss we shall as- 
sume to be at x Q = 0). The solution S(y) is obtained by 
iterative numerical algorithms, solvers, devised to invert 
so-called sparse matrices, like the matrices that result 
from the discretization of differential equations by finite 
differences methods. In this letter we shall not discuss 
the details of any particular solver (see ref. [T] for a com- 
plete review and for an updated list of references). For 
any solver one checks if the condition 



^2(D[U} + M) x , y S n (y)-rj(x) 



< r 



(2) 



is satisfyed. Here S n (y) is the tentative solution at itera- 
tion number n, the norm is any good norm in field space 
and r, the residue, is the global numerical accuracy re- 
quested for the solution. Typically r is a small number 
of the order of the arithmetic precision allowed by the 
computer architecture. Depending on the values of the 
quark mass the solution of eq. ([IJ poses different numeri- 
cal problems. For light quarks the matrix (£>[£/] + M) X:V 
is badly conditioned and its numerical inversion requires 
a big number of iterations. At the other extreme, the 
number of iterations required for heavy quark masses 
is small but there may be problems with the numerical 
accuracy resulting for the time-slices far away from the 
source (z/o 3> 0) . Indeed eq. ^ is a global condition while 
for heavy quark propagators the time-slices far away from 
the source are exponentially suppresed by a factor of the 
order of exp(— Myo) and give a negligible contribution to 
the norm on the left side of eq. pi). When this problem 



arises one cannot trust numerical results at large times 
and it becomes impossible to extract physical informa- 
tions by fitting the leading exponentials contributing to 
correlation functions. 

In order to alleviate both difficulties, we propose a pre- 
conditioning of the Dirac operator that factorises from 
the propagator a function aiming to modify its leading 
decay with the distance 1 . The simplest choice is to fac- 
torize a function a(yo), to solve numerically the precon- 
ditioned equation, and to restore the original propagator 
by multiplying each time slice for l/a(yo). a(yo) is de- 
fined such that all the different time-slices give compara- 
ble contributions to the calculation of the residue in the 
preconditioned case. Our preconditioning is inspired to 
what is usually done in deriving the Eichten and Hill [2 
lattice HQET action but of course does not introduce 
any approximation. Indeed, the choice above is suited for 
heavy quark propagators, while for light quark masses we 
will introduce a generalisation of the factorised function. 



II. PRECONDITIONING HEAVY QUARK 
PROPAGATORS 



We work with the (9(a)-improved Wilson lattice Dirac 
operator but the numerical problems that we address 
arise also with alternative discretisations of the contin- 
uum action and the proposed solution can as well be eas- 
ily implemented in those cases. We have tested our pre- 
conditioning scheme both for heavy and for light quark 
masses and in the free and in the interacting case. We 
start with the results for the heavy quarks. We first want 
to pick up a case where the problem arises. As an exam- 
ple, we have calculated the correlation function 



Cpp(yo) = [S\y)S(y)] 



(3) 



by solving eq. ([I]) for a heavy quark propagator of mass 
M ~ 0.5 in the free theory for different choices of the 



1 see refs. [3| [4] for different approaches 
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FIG. 1: The red points correspond to the correlation fuction — Cpp(yo) obtained by inverting the lattice Dirac equation in 
the free theory with M ~ 0.5 and with a residue r = 10" 11 . The black points correspond to the same quantity but have been 
obtained with a residue r = 10 -6 . The blue points correspond to the correlation function —C' PP (t/o) obtained by solving the 
preconditioned lattice Dirac equation with M ~ 0.5 and cto = 0.4. The two black lines correspond to r 2 for the two values of 
the residue used in the calculations. We use logarithmic scale on the y-axis. 
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FIG. 2: The red and black sets of points are the effective masses of the corresponding correlators shown in FIG. ??. The 
blue set of points is the effective mass of the corresponding correlator in the top panel multiplyied for the restoration factor as 
explained in the text. 



residue. More precisely the red points in FIG. [T] have 
been obtained with a residue r = 10 -11 while the black 
points with a residue r = 10~ 6 and the two black lines 
correspond to the squares of these two values of r. As 
is clearly visible from FIG. [TJ and from FIG. [2] where 
we show the effective masses of the correlations shown 
in FIG. [I] the black points start to deviate from the red 
ones for yo ~ 18, i.e. when the correlator, which in this 
case is just the square module of the propagator, becomes 
smaller than the square of the "loose" residue r = 10~ 6 . 

If the time extent of the lattice is not too large the 
problem can be solved by brute force by lowering the 
residue and the results obtained in the preset case with 
r = 10~ n can be considered as exact. If instead 
the time extent of the lattice is rather large the brute 
force approach cannot be considered because the required 



residues would be smaller than what is allowed on double- 
precision architectures, even in the case of moderately 
heavy quarks. In the case under consideration, by choos- 
ing a loose precision, i.e. a residue r = 10~ 6 , we make the 
numerical problem evident and we show that also such 
an "extreme" situation can be recovered by using our 
proposal. Moreover, we notice that a residue r = 10~ 6 
is the smallest allowed on single-precision architectures 
that presently are considerably much faster than double- 
precision ones. 

We now come to the proposed solution. We redefine 
the quark fields and the propagators as follows 

S(y) = a(y ) S'(y) 

v(y) = "(yo) v'(y) (4) 
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Once the previous expressions are inserted in eq. ([!]) we 
get the preconditioned system 



J2(D'[U]+M) x . y S'(y)=r 1 '(x) 



(5) 



that we solve numerically in place of eq. ([I]). In order to 
write the preconditioned Dirac operator it is sufficient 
to modify the forward and backward lattice covariant 
derivatives in the time direction accoring to 

V S(y) = U (y)S(y + 6)-S(y) 



_ °^+Au ( y )S'(y + 6)-S'(y) 



<x(yo) 



VlS(y) = S(y) 



S'(y) 



Ul{y 



6)S(y-6) 



<x(ya) 



Ut(y-0)S'(y~0) 



(6) 



Particular care has to be used at the boundaries of the 
lattice in order to respect the boundary conditions orig- 
inally satisfied by the quark fields. If as in the case of 
FIG. [I] S(y) satisfies anti-periodic boundary conditions 
along the time direction, it follows from eq. Q that 

S(y + L d) = -S(y) 



In FIG.[3]we show the same plot as in FIG.[2]but in the 
interacting theory. The gauge ensamble used correspond 
to the entry Eh in TABLE [I] The size of the lattice is 
L0L1L2L3 = 64 x 32 3 and the hopping parameter of the 
sea quarks is k sea — 0.13625 corresponding to a PCAC 
quark mass of about amf^ AC ~ 0.07. The data shown in 
FIG. [3] correspond to a pseudoscalar-pseudoscalar corre- 
lator, as in the free theory case, of two degenerate heavy 
quarks with hopping parameters kh = 0.125 correspond- 
ing to a PCAC quark mass of about am^ CAC ~ 0.35. 
The unpreconditioned correlators decay approximately 
as fast as in the free theory case and from the differ- 
ence of the black (unpreconditioned, r = 10~ 6 ) and red 
(unpreconditioned, r = 10~ ) sets of data we see the 
same distortion of FIG. [2] The blue points have been 
obtained by solving eq. (pi after having factorized a(yo) 
with ao = 0.4 and by restoring the results according to 
eq. ( 10 1 . Also in the interacting theory the blue points are 



identical to red points though they have been obtained 
with the same loose residue r = 10 -6 used to obtain the 
black points. 

We close this section by observing that our precondi- 
tioning technique may be particularly useful when work- 
ing with the Schrodinger Functional 5, 6 formulation of 
the theory. In this case, countrary to the case of periodic 
boundary conditions along the time direction, the cor- 
relators decay exponantially over the whole time extent 
of the lattice and one has to choose very small residues 
also in computing relatively light quark propagators. We 
have performed several succesful experiments with our 
preconditioning technique also in the Schrodinger Func- 
tional case by using a(xo) = exp(— ctQXo). 



S'(y + L b) 



a (2/0) 



a (2/0 + L 0) 



S'(y) 



(7) 



The blue points in FIG.[T]correspond to the correlation 
function 

^Pp(yo)--E tr [( 5 ') f ^) 5 '(y)] ( 8 ) 

v 

obtained by solving eq. ([5| with the loose residue r = 
10 -6 but after having factorized the function 

a(y ) = cosh[a (yo - W 2 )] ( 9 ) 

by setting «o = 0.4. As expected, the preconditioned cor- 
relator stays above the line of the loose precision residue 
and the "exact" result can be back recovered as follows 



Cpp(y ) = [a(yo)} 2 Cpp{y ) 



(10) 



In FIG.[2]the blue points correspond to the effective mass 
of the preconditioned correlator after the "restoration" of 
eq. ( 10 1 and fall exactly on top of the red ones in spite 



of the fact that they have been obtained with the same 
loose precision that affected the non preconditioned black 
points. 



III. PRECONDITIONING LIGHT QUARK 
PROPAGATORS 



In this section we shall briefly discuss how the ideas 
developed and discussed in the previous section can be 
used to accelerate the numerical calculation of light quark 
propagators. We start our discussion by generalizing 
eq. Q as follows 

S{y) = P(yo> 2/1.2/2,2/3) S'(y) 

v(v) = £(2/0,2/1,2/2,2/3) v'(v) (11) 

In the following we shall consider the particular choice 



/3(2/o,2/i,2/2,2/3) = n a{Vp) 



n 



£i coBh[a (^ - 



(12) 



and the preconditioned lattice Dirac operator can be ob- 
tained as easily as before by changing all the covariant 
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FIG. 3: The red points correspond to the effective mass of the correlation fuction —Cpp(yo) obtained by inverting the lattice 
Dirac equation in the interacting theory with ara' m ~ 0.35 and with a residue r = 10 -11 . The black points correspond to 
the same quantity but have been obtained with a residue r = 10~ 6 . The blue points correspond to the effective mass of the 
restored correlation function —C' PP (j/o) obtained by solving the preconditioned lattice Dirac equation with am^ ac ~ 0.35 and 
Qo = 0.4. 



j3 L0L1L2L3 k s 



«o iterations 
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TABLE I: Gauge configurations have been generated with 
nf — 2 dynamical 0(a)-improved Wilson quarks with c sw = 
1.90952. The figures in the last column correspond to the av- 
erage of the number of iterations required to invert the Dirac 
equation in the unitary theory by using the SAP+GCR in- 
verter for several values of the preconditioning parameter ao • 
The values corresponding to ao — 0.0 have been obtained 
without using our preconditioning technique. 



derivatives according to 

V^(y) - U^y)S(y + fi)-S(y) 

a(y M + 1) 



U^y)S'(y + fl)-S'(y) 



VtS(y) = S(y)-Ul(y-p,)S(y-fi) 



H. S'(y)- a{y » ^ UKy-ttS'iy-d) 



(13) 



and by changing accordingly the boundary conditions in 
all directions as done in eqs. Q for the time direction. 

An important difference of the present case with re- 
spect to the one discussed in the previous section is that 
the restoration of the true propagator must be performed 
before making the contractions needed to build correla- 
tion functions by using the first of eqs. (111. 



Here the preconditioning is not to gain precision, but 
to accelerate the convergence of the inversion. Therefore, 
by applying eqs. ( 11 ), ( 12 1 and ( 13 ) to the calculation of a 



light quark propagator one aims to make the propagator 
to decay faster than the original unpreconditioned oper- 
ator. By judiciously chosing the parameter ao it is possi- 
ble to change the condition number of the preconditioned 
system without compromising the numerical accuracy of 
the solution, an operation that should be performed on 
double-precision computer architectures. 

In TABLE [T] we quantify the gain in computational 
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time that can be achieved by showing the number of it- 
erations of the SAP+GCR solver required to solve the 
lattice Dirac equation for light quarks with and without 
our preconditioning. The SAP+GCR solver has been in- 
troduced and explained in details by the author in ref . [7] . 
The gauge ensambles used to perform this test have been 
generated within the CLS agreement [H] with the param- 
eters given in the table. In the case of the _E-lattices the 
SAP+GCR solver has been ran on 128 processors of a 
cluster of PC's by dividing the global lattices into blocks 
of 4 4 points. In the case of the D-lattice the SAP+GCR 
solver has been ran on 32 processors of a cluster of PC's 
by dividing the global lattices into blocks of 6 4 points. 
The table shows that by increasing the value of the pa- 
rameter a the number of iterations goes down with a 
time gain that can easily reach the 30%. In the case 



under discussion, we checked that higher values of ao 
would induce a "heavy quark" like behavior and produce 
distorted results for the reasons discussed at lenghty in 
the previous section. 

The method discussed in this letter can be generalised 
by adding some Dirac structure in the factorised function, 
an option presently under investigation. 
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